load("lupus_data_analysis.RData")
#library(SparseM)
source("smoother.r")

####### Read in Lupus Data ##############
load("steve.RData")
#steve_fan<- final_all

toronto<- final_all[which(!is.na(final_all$CTUID)),]
rm(list=c("steve_fan", "final_all"))
toronto$Observed<- replace(toronto$Observed, is.na(toronto$Observed), 0)

pixel<- unique(toronto[c(2, 3)])
pixel<- pixel[order(pixel$x, pixel$y),] 
pixel<- data.frame(pixel, PID=1:nrow(pixel), centroid=complex(real=pixel$x, imaginary=pixel$y))

region<- unique(toronto[c(4, 5)])
region<- region[order(region$Year, region$CTUID), ]
region<- data.frame(region, RID=1:nrow(region))

toronto<- merge(toronto, pixel)
toronto<- merge(toronto, region)
toronto<- toronto[order(toronto$x, toronto$y),]

attach(toronto)
incident.mat<- as.matrix.csr(0, nrow=nrow(region), ncol=nrow(pixel))
for (i in 1:nrow(region)) incident.mat[i, PID[which(RID==i)]]<- 1

tol<- 1E-8
cellwidth<- 175
#evalpt<- as.vector(outer(as.vector(outer(quadPts*(cellwidth/2), quadPts*(cellwidth/2), function(x, y) complex(real=x, imaginary=y))), centroid, "+"))

#######  construct the smoothing matrix ###############
#bwSeq<- seq(0, 200, by=5)
bwSeq<- c(0, seq(150, 1500, by=50))
offsets<- as.vector(tapply(Expected, RID, unique) %*% incident.mat)
offsets.mat<- as(offsets,"matrix.diag.csr")

for (j in (1:length(bwSeq))) {
  K.mat<- t(ems_constructor(pixel$centroid, offsets, bwSeq[j], cellwidth, quadPts, quadWgts))
  save(K.mat, file=paste("kernel.matrix", j, "RData", sep="."))
  cat(j, "\n")
}


#kernel.fun<- function(x) {
#  x<- complex(real=replace(Re(x), is.nan(Re(x)), 0), imaginary=replace(Im(x), is.nan(Im(x)), 0))
#  dist_sq<- Mod(x)^2;
#  output<- (3/pi)*(1-dist_sq)^2;
#  replace(output, dist_sq>1, 0);
#}
############ Observed Counts ####################
observed<- tapply(Observed, RID, unique)
expected<- tapply(toronto$Expected, toronto$RID, unique)

#################  EMS #########################################
attach(toronto)
niter.mat<- rep(0, length(bwSeq))

output<- pixel
for (j in (1:length(bwSeq))) output<- cbind(output, est=rep(0, nrow(pixel)))
names(output)[5:33]<- paste("est", 1:length(bwSeq), sep=".")

for (j in (2:length(bwSeq))) {
#for (j in (28:29)) {
	load(file=paste("kernel.matrix", j, "RData", sep="."))
	##################
	# the diagnoal matrix whose diagnoal entry is Lambda.old
	Lambda.old<- as.matrix.csr(matrix(1, nrow=1, ncol=length(pixel$centroid)))
	tmp.old.mat<- as(as.vector(as.matrix(Lambda.old)), "matrix.diag.csr")
	####################
	niter<- as.integer(1)
      
	repeat {
		##### EM step #####
		V.inv<- diag(1/as.vector(as.matrix(Lambda.old %*% t(incident.mat))))
		dd<- observed %*% V.inv
		dd<- replace(dd, is.nan(dd), 0)
		tmpEM<- dd %*% incident.mat %*% tmp.old.mat %*% offsets.mat
		#tmpEM<- (observed %*% V.inv) %*% incident.mat %*% tmp.old.mat %*% offsets.mat
		##### S step #####
		Lambda<-  tmpEM %*% (offsets.mat %*% K.mat)
		check<- mean((as.vector(as.matrix(Lambda - Lambda.old)))^2)
		if (check > tol) {
			cat(niter, "\t", check, "\n")
			Lambda.old<- Lambda
			##################
			# the diagnoal matrix whose diagnoal entry is Lambda.old
			tmp.old.mat<- as(as.vector(as.matrix(Lambda.old)),"matrix.diag.csr")
			####################
			niter<- as.integer(1+niter);
		}
		else {
			niter.mat[j]<- niter
			output[,(4+j)]<- as.vector(as.matrix(Lambda))
			break
		}
	}
}
result<- merge(pixel, output)
rm(list=c("dd", "tmp.old.mat", "check", "tmpEM", "V.inv", "niter", "K.mat"))

pdf(file="lupus_map_with_exact_obs_oct282008.pdf")
levelplot(est.1 ~ x*y, data=output, cuts=20, region=TRUE, col.regions=gray(1:50/50)[21:40], colorkey=list(col=gray(1:50/50)[21:40]))
trellis.focus(); ltext(exact$x, exact$y, "x")
levelplot(est.1 ~ x*y, data=output, cuts=20, region=TRUE, col.regions=gray(1:50/50)[21:40], colorkey=list(col=gray(1:50/50)[20:40]), xlim=c(620000, 640000), ylim=c(4825000, 4845000), main="EM Relative Risk Estimate")
trellis.focus(); ltext(exact$x, exact$y, "x")
levelplot(est.5 ~ x*y, data=output, cuts=20, region=TRUE, col.regions=gray(1:50/50)[21:40], colorkey=list(col=gray(1:50/50)[20:40]), xlim=c(620000, 640000), ylim=c(4825000, 4845000), main=paste("EMS Relative Risk Estimate with the Bandwidth of", bwSeq[5]))
trellis.focus(); ltext(exact$x, exact$y, "x")
levelplot(est.10 ~ x*y, data=output, cuts=20, region=TRUE, col.regions=gray(1:50/50)[21:40], colorkey=list(col=gray(1:50/50)[20:40]), xlim=c(620000, 640000), ylim=c(4825000, 4845000), main=paste("EMS Relative Risk Estimate with the Bandwidth of", bwSeq[10]))
trellis.focus(); ltext(exact$x, exact$y, "x")
levelplot(est.15 ~ x*y, data=output, cuts=20, region=TRUE, col.regions=gray(1:50/50)[21:40], colorkey=list(col=gray(1:50/50)[20:40]), xlim=c(620000, 640000), ylim=c(4825000, 4845000), main=paste("EMS Relative Risk Estimate with the Bandwidth of", bwSeq[15]))
trellis.focus(); ltext(exact$x, exact$y, "x")
levelplot(est.20 ~ x*y, data=output, cuts=20, region=TRUE, col.regions=gray(1:50/50)[21:40], colorkey=list(col=gray(1:50/50)[20:40]), xlim=c(620000, 640000), ylim=c(4825000, 4845000), main=paste("EMS Relative Risk Estimate with the Bandwidth of", bwSeq[20]))
trellis.focus(); ltext(exact$x, exact$y, "x")
levelplot(est.25 ~ x*y, data=output, cuts=20, region=TRUE, col.regions=gray(1:50/50)[21:40], colorkey=list(col=gray(1:50/50)[20:40]), xlim=c(620000, 640000), ylim=c(4825000, 4845000), main=paste("EMS Relative Risk Estimate with the Bandwidth of", bwSeq[25]))
trellis.focus(); ltext(exact$x, exact$y, "x")
levelplot(est.26 ~ x*y, data=output, cuts=20, region=TRUE, col.regions=gray(1:50/50)[21:40], colorkey=list(col=gray(1:50/50)[20:40]), xlim=c(620000, 640000), ylim=c(4825000, 4845000), main=paste("EMS Relative Risk Estimate with the Bandwidth of", bwSeq[26]))
trellis.focus(); ltext(exact$x, exact$y, "x")
dev.off()

pdf(file="lupus_map_oct282008.pdf")
levelplot(est.1 ~ x*y, data=output, cuts=20, region=TRUE, col.regions=gray(1:50/50)[21:40], colorkey=list(col=gray(1:50/50)[20:40]), xlim=c(620000, 640000), ylim=c(4825000, 4845000), main="EM Relative Risk Estimate")
levelplot(est.5 ~ x*y, data=output, cuts=20, region=TRUE, col.regions=gray(1:50/50)[21:40], colorkey=list(col=gray(1:50/50)[20:40]), xlim=c(620000, 640000), ylim=c(4825000, 4845000), main=paste("EMS Relative Risk Estimate with the Bandwidth of", bwSeq[5]))
levelplot(est.10 ~ x*y, data=output, cuts=20, region=TRUE, col.regions=gray(1:50/50)[21:40], colorkey=list(col=gray(1:50/50)[20:40]), xlim=c(620000, 640000), ylim=c(4825000, 4845000), main=paste("EMS Relative Risk Estimate with the Bandwidth of", bwSeq[10]))
levelplot(est.15 ~ x*y, data=output, cuts=20, region=TRUE, col.regions=gray(1:50/50)[21:40], colorkey=list(col=gray(1:50/50)[20:40]), xlim=c(620000, 640000), ylim=c(4825000, 4845000), main=paste("EMS Relative Risk Estimate with the Bandwidth of", bwSeq[15]))
levelplot(est.20 ~ x*y, data=output, cuts=20, region=TRUE, col.regions=gray(1:50/50)[21:40], colorkey=list(col=gray(1:50/50)[20:40]), xlim=c(620000, 640000), ylim=c(4825000, 4845000), main=paste("EMS Relative Risk Estimate with the Bandwidth of", bwSeq[20]))
levelplot(est.25 ~ x*y, data=output, cuts=20, region=TRUE, col.regions=gray(1:50/50)[21:40], colorkey=list(col=gray(1:50/50)[20:40]), xlim=c(620000, 640000), ylim=c(4825000, 4845000), main=paste("EMS Relative Risk Estimate with the Bandwidth of", bwSeq[25]))
dev.off()


### cross validation ###
#region.pixel<- unique(toronto[c(1, 8:10)])
#xyplot(Im(centroid) ~ Re(centroid)|Year, data=region.pixel)

tmp.offsets<- tapply(toronto$Expected, toronto$RID, unique)
tmp.RID<- tapply(toronto$RID, toronto$Year, function(x) unique(sort(x)))
#tmp.RID[[4]]<- union(which(tmp.offsets==0), tmp.RID[[4]])

pixel.ex<- list()
incident.mat.cv<- list()
offsets.cv.lst<- offsets.cv.mat.lst<- list()
for (i in 1:4) {
  ############ pixel excluded if the ith map is left out ####################
  pixel.ex[[i]]<- which(0==as.vector(tmp.offsets[-tmp.RID[[i]]] %*% incident.mat[-tmp.RID[[i]],]))
  ############# The incident matrix for cross validatioin ####################
  incident.mat.cv[[i]]<- incident.mat[-tmp.RID[[i]], -pixel.ex[[i]]]
  ############# Offsets ######################################################
  offsets.cv.lst[[i]]<- as.vector(tmp.offsets[-tmp.RID[[i]]] %*% incident.mat.cv[[i]])
  offsets.cv.mat.lst[[i]]<- as(offsets.cv.lst[[i]], "matrix.diag.csr")
}

save.image("lupus_data_analysis_new_data.RData")


load("lupus_data_analysis_new_data.RData")
library(SparseM)
library(lattice)

############## Smoothing Matrices ###########
bwSeq.cv<- 150*(0:15)

#for (j in (1:length(bwSeq.cv))) {
#for (j in (1:7)) {
for (j in 1:15) {
  K.mat.cv<- list()
  K.mat.cv[[1]]<- t(ems_constructor(pixel$centroid[-pixel.ex[[1]]], offsets.cv.lst[[1]], bwSeq.cv[j], cellwidth, quadPts, quadWgts)); gc()
  K.mat.cv[[2]]<- t(ems_constructor(pixel$centroid[-pixel.ex[[2]]], offsets.cv.lst[[2]], bwSeq.cv[j], cellwidth, quadPts, quadWgts)); gc()
  K.mat.cv[[3]]<- t(ems_constructor(pixel$centroid[-pixel.ex[[3]]], offsets.cv.lst[[3]], bwSeq.cv[j], cellwidth, quadPts, quadWgts)); gc()
  K.mat.cv[[4]]<- t(ems_constructor(pixel$centroid[-pixel.ex[[4]]], offsets.cv.lst[[4]], bwSeq.cv[j], cellwidth, quadPts, quadWgts)); gc()
  save(K.mat.cv, file=paste("kernel.matrix.cv", j, "RData", sep="."))
  cat(j, "\n")
}

region.pixel<- unique(toronto[c(1, 8:10)])
output.cv<- output.cv.lst<- list()

for (j in (1:length(bwSeq))) {
  tmp.output<- pixel[c(4, 3)]
  load(file=paste("kernel.matrix.cv", j, "RData", sep="."))
  for (i in 1:4) {
    ##################
    # the diagnoal matrix whose diagnoal entry is Lambda.old
    #	Lambda.old<- observed %*% incident.mat %*% as(1/as.vector(rep(1, nrow(region)) %*% incident.mat), "matrix.diag.csr")
    Lambda.old<- as.matrix.csr(1, nrow=1, ncol=length(pixel$centroid[-pixel.ex[[i]]]))
    tmp.old.mat<- as(as.vector(as.matrix(Lambda.old)), "matrix.diag.csr")
    ####################
    repeat {
	##### EM step #####
	V.inv<- diag(1/as.vector(as.matrix(Lambda.old %*% t(incident.mat.cv[[i]]))))
	dd<- observed[-tmp.RID[[i]]] %*% V.inv
	dd<- replace(dd, is.nan(dd), 0)
	tmpEM<- dd %*% incident.mat.cv[[i]] %*% tmp.old.mat %*% offsets.cv.mat.lst[[i]]
	gc()
	##### S step #####
	Lambda<-  tmpEM %*% (offsets.cv.mat.lst[[i]] %*% K.mat.cv[[i]])
	#Lambda<-  tmpEM %*% (offsets.cv.mat.lst[[i]] %*% K.mat.cv)
	gc()
	check<- mean((as.vector(as.matrix(Lambda - Lambda.old)))^2)
	if (check > tol) {
		cat(check, "\n")
		Lambda.old<- Lambda
		##################
		# the diagnoal matrix whose diagnoal entry is Lambda.old
		tmp.old.mat<- as(as.vector(as.matrix(Lambda.old)), "matrix.diag.csr")
		####################
	}
	else {
	  Lambda<- as.vector(as.matrix(Lambda))
	  dd<- data.frame(PID=tmp.output$PID[-pixel.ex[[i]]], est=Lambda)
	  output.cv[[i]]<- tmp.output<- merge(tmp.output, dd, all.x=TRUE)
	  break
	}
      }
    #output.cv[[i]]<- tmp.output
    }

    output.cv.lst[[j]]<- output.cv
}

rm(list=c("dd", "tmp.old.mat", "check", "tmpEM", "V.inv", "tmp.output", "K.mat.cv"))




pred.error<- rep(0, 15)
for (j in 1:15) {
  for (i in 1:4) {
    #tmp<- as(1/as.vector(rep(1, nrow(pixel)) %*% t(incident.mat[tmp.RID[[i]],])), "matrix.diag.csr") %*% incident.mat[tmp.RID[[i]],]
    #intRR<- as.vector(replace(output.cv.lst[[j]][[i]]$est, is.na(output.cv.lst[[j]][[i]]$est), 0) %*% t(tmp) )
    #pred.error[j]<- pred.error[j] + mean((observed[tmp.RID[[i]]] - expected[tmp.RID[[i]]]*intRR)^2)
    intRR<- replace(output.cv.lst[[j]][[i]]$est, is.na(output.cv.lst[[j]][[i]]$est), 0); #intRR<- intRR/max(intRR)
    intRR<- as.vector( intRR %*% t(incident.mat[tmp.RID[[i]],]) )
    pred.error[j]<- pred.error[j] + sum((observed[tmp.RID[[i]]] - expected[tmp.RID[[i]]]*intRR)^2)/4
  }
  print(pred.error)
}

bwSeq.cv[which(min(pred.error)==pred.error)] 
# 750


which(bwSeq==1350) # j=26

pdf(file="lupus_map_using_CV.pdf")
levelplot(est.1 ~ x*y, data=output, cuts=20, region=TRUE, col.regions=gray(1:50/50)[21:40], colorkey=list(col=gray(1:50/50)[21:40]))
trellis.focus(); ltext(exact$x, exact$y, "x")
levelplot(est.1 ~ x*y, data=output, cuts=20, region=TRUE, col.regions=gray(1:50/50)[21:40], colorkey=list(col=gray(1:50/50)[20:40]), xlim=c(620000, 640000), ylim=c(4825000, 4845000), main="EM Relative Risk Estimate", sub="Area-censored Data")
trellis.focus(); ltext(exact$x, exact$y, "x")
levelplot(est.26 ~ x*y, data=output, cuts=20, region=TRUE, col.regions=gray(1:50/50)[21:40], colorkey=list(col=gray(1:50/50)[20:40]), xlim=c(620000, 640000), ylim=c(4825000, 4845000), main=paste("EMS Relative Risk Estimate with the Bandwidth of", bwSeq[26]), sub="Area-censored Data")
trellis.focus(); ltext(exact$x, exact$y, "x")



levelplot(est.26 ~ x*y, data=output.exact[which(output.exact$x>=620000 & output.exact$x<=640000 & output.exact$y>=4825000 & output.exact$y<=4845000), ], cuts=20, region=TRUE, col.regions=gray(1:50/50)[21:40], colorkey=list(col=gray(1:50/50)[20:40]), xlim=c(620000, 640000), ylim=c(4825000, 4845000), main=paste("Relative Risk Estimate with the Bandwidth of", bwSeq[26]), sub="Exact Data")
trellis.focus(); ltext(exact$x, exact$y, "x")
levelplot(est.26 ~ x*y, data=output.exact, cuts=20, region=TRUE, col.regions=gray(1:50/50)[21:40], colorkey=list(col=gray(1:50/50)[20:40]), main=paste("Relative Risk Estimate with the Bandwidth of", bwSeq[26]), sub="Exact Data")
trellis.focus(); ltext(exact$x, exact$y, "x")
#levelplot(est.1 ~ x*y, data=output, cuts=20, region=TRUE, col.regions=gray(1:50/50)[21:40], colorkey=list(col=gray(1:50/50)[20:40]), xlim=c(620000, 640000), ylim=c(4825000, 4845000), main="EM Relative Risk Estimate")
#levelplot(est.26 ~ x*y, data=output, cuts=20, region=TRUE, col.regions=gray(1:50/50)[21:40], colorkey=list(col=gray(1:50/50)[20:40]), xlim=c(620000, 640000), ylim=c(4825000, 4845000), main=paste("EMS Relative Risk Estimate with the Bandwidth of", bwSeq[26]))
dev.off()





commonPixel<- pixel[which(4==as.vector(rep(1, 2783) %*% incident.mat)),]
commonPID<- which(4==as.vector(rep(1, 2783) %*% incident.mat))
dd<- cbind(as.vector(rep(1, 206957)%*%t(incident.mat)), as.vector(rep(1, 100930) %*% t(incident.mat[,commonPID])))
commonRID<- which(dd[,1]>0 & dd[,2]>0)


pred.error.common<- rep(0, 16)
for (j in 1:16) {
  for (i in 1:4) {
    intRR<- replace(output.cv.lst[[j]][[i]]$est, is.na(output.cv.lst[[j]][[i]]$est), 0)
    intRR<- as.vector( intRR %*% t(incident.mat[intersect(tmp.RID[[i]], commonRID),]) )
    pred.error.common[j]<- sum((observed[intersect(tmp.RID[[i]], commonRID)] - expected[intersect(tmp.RID[[i]], commonRID)]*intRR)^2)/4 + pred.error.common[j]
  }
  print(pred.error.common)
}








### exact ###
#load("exact.RData")
output.exact<- pixel
for (j in (1:length(bwSeq))) output.exact<- cbind(output.exact, est=rep(0, nrow(pixel)))
names(output.exact)[5:33]<- paste("est", 1:length(bwSeq), sep=".")

location<- complex(real=exact$x, imaginary=exact$y)
for (j in 1:length(bwSeq)) {
  tmpD<- intEstimate.exact(location, pixel$centroid, offsets, bandwidth=bwSeq[j], cellwidth, quadPts, quadWgts)
  output.exact[, (4+j)]<- tmpD$Lambda
  cat(j, "\n")
}
########### cross validation #####


bwSeq.cv<- 150*(7:20)
output.cv.exact.lst<- output.cv.exact<- list()
census.year<- unique(exact$Year)
for (j in 12:length(bwSeq.cv)) {
  for (i in 1:4) {
    tmplocation<- complex(real=exact$x[which(exact$Year!=census.year[i])], imaginary=exact$y[which(exact$Year!=census.year[i])])
    tmpD<- intEstimate.exact(tmplocation, pixel$centroid, offsets, bandwidth=bwSeq.cv[j], cellwidth, quadPts, quadWgts)
    output.cv.exact[[i]]<- tmpD$Lambda 
  }
  output.cv.exact.lst[[j]]<- output.cv.exact
  cat(j, "\n")
  rm(list=c("tmplocation", "tmpD"))
  save.image("lupus_data_analysis_new_data.RData")
}

rm(list=c("tmplocation", "tmpD"))
save.image("lupus_data_analysis_new_data.RData")


## give the pixel that contain the disease incident ####
dd<- outer(location, pixel$centroid, function(x, y) Mod(x-y))
evtPID<- apply(dd, 1, function(x) which(x==min(x)))

## to cross validate, we need to get the expected number of events for the pixels that contain events in a census period ##
tmpPixel<- merge(pixel, unique(toronto[c(1, 6, 8, 9)]))
tmpPixel<- tmpPixel[order(tmpPixel$Year, tmpPixel$PID),]


pred.error.exact<- rep(0, length(bwSeq.cv))
for (j in 1:6) {
  for (i in 1:4) {
    exp.spatial<- output.cv.exact.lst[[j]][[i]][sort(unique(evtPID[which(exact$Year==census.year[i])]))]/(cellwidth^2)
    exp.cov<- tmpPixel$Expected[which(tmpPixel$Year==census.year[i] & tmpPixel$PID %in% evtPID[which(exact$Year==census.year[i])])]
    pred.error.exact[j]<- pred.error.exact[j] + sum((table(evtPID[which(exact$Year==census.year[i])]) - exp.spatial*exp.cov)^2)/4
  }
}

rm(list=c("tmpPixel", "exp.spatial", "exp.cov", "evtPID", "dd"))
save.image("lupus_data_analysis_new_data.RData")